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Abstract. We present a detailed mathematical analysis of the Fourier response of binary pulsar signals whose 
frequencies are modulated by circular orbital motion. The fluctuation power spectrum of such signals is found to 
be i^orb-periodic over a compact frequency range, where forb denotes orbital frequency. Subsequently, we consider 
a wide range of binary systems with circular orbits and short orbital periods, and present a Partial Coherence 
Recovery Technique for searching for binary millisecond X-ray and radio pulsars. We use numerical simulations to 
investigate the detectability of pulsars in such systems with Porb< 6 hours, using this technique and three widely 
used pulsar search methods. These simulations demonstrate that the Partial Coherence Recovery Technique is 
on average several times more sensitive at detecting pulsars in close binary systems when the data span is more 
than 2 orbital periods. The systems one may find using such a method can be used to improve the constraints 
on the coalescence rate of compact objects and they also represent those systems most likely to be detected with 
gravitational wave detectors such as LISA. 
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1. Introduction 

The importance of binarity and the connection with the 
low-mass X-ray binaries (LMXBs) in the formation of 
millisecond pulsars (Alpar et al. 1982; Radhakrishnan & 
Srinivasan 1982) extends back to the discovery of the first 
millisecond pulsar (Backer et al. 1982). In the intervening 
years, a large number of binary millisecond radio pulsars 
and X-ray pulsars have been discovered (for recent reviews 
see Kramer et al. 2000). However, these tend to be either 
bright pulsars, members of wide binaries or to have low- 
mass companions. The remaining systems are more diffi- 
cult to detect owing to large orbital accelerations which 
cause, via the time-dependent Doppler effect, the apparent 
pulse period to change significantly during long integra- 
tions. Compensating completely for this smearing of the 
pulsed signal is extremely computationally expensive as a 
large range of orbital parameters has to be searched. 
Numerous authors have developed and used techniques 
to partially correct for this smearing, e.g. Middleditch & 
Priedhorsky (1986); Anderson et al. (1990); Wood et al. 
(1991); Ransom (1999). In 1990 Anderson etai. Anderson 
et al. (1990) used a constant acceleration technique to suc- 
cessfully discover PSR 2127-l-llC in the globular cluster 
M15. Despite the success of this survey, it was not un- 
til recently that the true extent of the binary millisecond 
pulsar population that could be uncovered by using such 
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techniques became apparent (Camilo et al. 2000; D'Amico 
et al. 2000; Ransom et al. 2000). However a whole possi- 
ble population of rapidly rotating radio/X-ray pulsars in 
close binary systems still remains to be investigated as 
even the current best search strategies are not very sen- 
sitive to them. This population includes faint pulsars in 
similar type binaries, pulsars in tighter binaries and those 
with more massive companions. 

Following the discovery of a large number of binary mil- 
lisecond pulsars in the globular cluster 47Tuc (Camilo 
et al. 2000), Rasio, Pfahl & Rappaport (2000) tried to 
model this observed population. They find that a large 
number of neutron star - white dwarf binaries in even 
tighter orbits and possibly more massive companions than 
those already discovered should exist in 47Tuc. The ex- 
istence of very short orbital period LMXBs such as 4U 
1820-30 (Stella, White & Priedhorsky 1987), and mil- 
lisecond pulsars with relatively massive companions, e.g. 
PSR B1744-24A (Lyne et al. 1990), also indicate that 
very short period binary millisecond pulsars should ex- 
ist in globular clusters (Bisnovatyi-Kogan 1989; Ergma & 
Fedarova 1991). Furthermore, multi-frequency radio imag- 
ing of globular clusters has revealed several unidentified 
steep-spectrum radio sources whose variable flux-densities 
and spectral indexes suggest radio pulsars (Fruchter & 
Goss 1990; 2000). Subsequent pulsed searches (e.g. Lyne 
et al. 1990; D'Amico 1993; Biggs et al. 1994) have de- 
tected millisecond pulsars coincident with some of these 
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point sources supporting the idea that the remainder are 
also pulsars which are so far undetected. Moreover, the 
number of pulsars detected in globular clusters is far less 
than the expected population based on the integrated flux 
density of the cores of some of the clusters. It has been 
argued that even bright pulsars could have been missed 
by previous searches due to Doppler smearing caused by 
the orbital motion in close binary systems, even though 
their radio emission remains detectable by interferomet- 
ric observations. A further source of possible pulsar can- 
didates are the subset of the EGRET unidentified 7-ray 
sources whose distribution above the Galactic disk has a 
scale height of about 2 kpc (Grenier 2001). These sources 
have so far evaded detection as radio pulsars and this may 
be because they are members of close binaries. 
The physics of gravitation in such extreme systems has re- 
cently received growing attention because such binary sys- 
tems will be the most common known continuous sources 
of gravitational radiation for gravitational wave detectors. 
The Laser Interferometer Space Antenna (LISA) will be 
especially sensitive to compact neutron star binaries with 
ultra-short orbital periods and large (> 0.5 Mq) compan- 
ion masses (Benacquista, Portegies Zwart & Rasio 2000). 
Detecting and studying pulsed radio or X-ray emission 
from the neutron stars in these systems will provide bet- 
ter constraints on the binary parameters and therefore 
be of paramount importance in helping to understand the 
gravitational wave emission (Dhurandhar & vecchio 2000) . 
Discovering such systems is further motivated by improv- 
ing the constraints on the coalescence rate of compact 
objects which remains difficult to estimate. Therefore, 
both theoretical predictions and observational evidence 
strongly encourage searching for highly accelerated mil- 
lisecond pulsars in globular clusters. 

Thus motivated, we present here a detailed mathemati- 
cal analysis of the Fourier response of a radio/X-ray bi- 
nary pulsar signal. We then discuss the sensitivity of three 
widely used search techniques applied to a range of binary 
systems. In Sec. ^, we present a complete description of 
our Partial Coherence Recovery Technique, which is based 
on the Phase Modulation Searching method first discussed 
by Ransom (Ransom 1999). We also discuss the computa- 
tional cost of the method, how to estimate detection lev- 
els and how to derive all the binary parameters directly 
from the Fourier signature of the binary. Finally we com- 
pare the sensitivity of these different search techniques 
and show that our method is not only computationally 
very efficient but also significantly more sensitive to very 
tight binaries. 

A list of the main symbols we use in this paper accompa- 
nied by a short description can be found in Appendix |^. 

2. Pulsar signal 

Consider a pulsar with a pulse frequency Vpsr, whose dis- 
tance from the observer d{t) is changing with time t. The 
pulsar signal s{t) emitted at time t is received by the ob- 



server at time t + d{t)/c, and can be represented by its 
harmonic decomposition 



5W = E*»W ' 5«(t)=A.(i)e'"^"''-'[*+'*(*)/^l (1) 

n=0 

where An{t) is the complex amplitude of the n*^ harmonic, 
c the velocity of light and j ~ We consider binary 

systems in circular orbits, for which the projected distance 
to the observer is 



d{t) — [ai sin(i)] cos (27ri^orbi + 0orb) + Constant (2) 

where ai is the pulsar's orbital radius, i the orbital incli- 
nation in radians, Voib the orbital frequency and 0orb the 
initial orbital phase measured from superior conjunction 
of the pulsar. For simplicity, we shall take Constant = 0. 
Expanding Eq. (§) in Eq. (|l|) yields 

where ipn = 2T:ni'ps,- [ai sin(i)/c] represents half the phase 
rotation in radians experienced by the n'^ harmonic of the 
pulsar signal during the orbit. Since the signal s„(t) is the 
n^^ harmonic of the Fourier series decomposition, we shall 
study its Fourier response. 

3. Fourier analysis 

We denote the Fourier response of any function ^ by \P 
and define it by ^{v) — J-[i>{t)]{i') where i/ indicates the 
frequency and J-' the Fourier Transform operator. As re- 
gards notations, Int[] shall define the integer part of and 
Frac[] the fractional part of. Henceforth for consistency 
with future notations, we call T-space what is usually re- 
ferred to as the time domain. 

In an ideal case, the complex amplitude An{t) = An, and 
the pulsar signal is a periodic function. However, in most 
astrophysical cases the periodicity of the signal is violated 
because 

1. The pulsed signal strength might vary during the ob- 
servation due to eclipses, interstellar scintillation and 
observational windowing. 

2. The pulsed signal could be modulated by stochastic 
or deterministic processes such as radio pulsar nulling, 
drifting subpulses or microstructure. 

The latter effects are not so important because the 
average pulse profile of radio pulsars is known to remain 
stable within the observing window. Nevertheless in prac- 
tice, An(t) is a function of time but the modulation itself 
remains to a large extent a priori unknown to the observer. 
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Fig. 1. (a) left: Jk{fn) with = 20. We plot fc e M whereas in reality /c e N. The Fourier response slightly extends 
beyond the amplitude of the phase modulation as e{ipn) 7^ 0. (b) right: in the general case, we plot s{(pn) as a function 

of (fn- 



3.1. Continuous Fourier response 



Let ^n{t) = 



s(2Triyorbt+0orb) define the phase rotation 



factor of Eq. and W{t) a real window function of fi- 
nite duration. For an observation constrained by a window 
W{t), Eq. (H) can be rewritten as 



(4) 



According to Fourier theory and by virtue of the shift 
theorem, 



(5) 



where * denotes the convolution operator. As we aim to 
study the Fourier response of the pulsating signal, we shall 
start with the calculation of ^n(i^). By definition. 



Ut)e~^^'''' dt. 



(6) 



The integration of Eq. (^) requires the substitution of ^n{t) 
by its series expansion given by (see Appendix p) 



(7) 



fc=-c 



where Jk are Bessel Functions of the First Kind of integer 
order and fc G N (see also Middleditch 1981). Combining 
Eqs. (|) and (g) yields 



OO 



jfc[0o,b+Tr/2] 



(5(j/ - fcl/orb) (8) 



where S denotes Dirac's delta function. This convolution 
greatly simplifies Eq. (N) into 



6.(^) = 







otherwise 



(9) 



In the ideal case of an infinite and continuous observa- 
tion, ^„ (ly) represents the Fourier response of the n*^ har- 
monic of a signal of constant amplitude emitted by a pul- 
sar whose phase is modulated by circular orbital motion. 
Following Eq. (||), ^„(i^) is defined in Eq. (||) relative to 
the n"^ harmonic frequency taken as the origin of the fre- 
quency scale. Of particular interest is its power spectrum, 
defined by P(^,n{v) = |^„(i/)p. A number of important 
properties of are : 

1. Periodicity 

P(^^n is a periodic sequence with period Vorh whose am- 
phtude is modulated by Jl{(pn). 

2. Compactness 

Since k is an integer, J-k{(pn) = (— l)'^ Jfc(¥'n) hence 
P^.n is symmetric with respect to fc = 0. Furthermore, 

3e{ipn) , y k > [ipn + e{ipn)] : 



\JkiVn)\ < 



Max[Jfe(^„)]fc 



1000 



(see Abramowitz & Stegun 1974 and Fig. |l|(a)). 
As shown in Fig. 0(b), e{(pn) < so Jk{fn) 
only slightly extends beyond its argument. Therefore, 
P^^n is compact on the interval /((/?«, forb) = 
[-E((/?„)i/orb;E((^„)i'orb] where E((^„) = ipn + e{ipn) 
radians. Note that 



k— — co 



fn t'orb 



(10) 



4 



Jouteux et al.: Searching for pulsars in close binary systems 




[Vorb/c] Npp^^ (Bins) 



Fig. 2. Top : pulsar radial velocity in units of Fourier fre- 
quency bins (vorb denotes the radial velocity of the pul- 
sar) . Bottom : power spectrum of the l*** harmonic shifted 
back to the origin of the frequency scale in arbitrary 

units. Here, ipi = 100, Np^^^^ = 2.0 and (port = 0. The re- 
sponse is symmetric because the pulsar signal is periodic 
in the observing window and -/Vp„,i, is an integer. 

where t^max denotes the maximum radial velocity of 
the pulsar. Thus, (/?„ i^orb represents the amplitude of 
the Doppler shift experienced by the n^^ harmonic of 
the pulsar. P^.„ is a comb of ~ 2(/?„ sidebands with a 
spacing of I'orb- 



3.2. Discrete Fourier response 

By necessity, any observed signal is of finite extent. The 
extent may be adjustable and selectable, but must be 
finite. Processing a finite-duration observation by esti- 
mating its complex spectrum directly from the Fourier 
Transform encounters many difficulties which make a bi- 
nary pulsar's Fourier response more complicated. 



3.2.1. Frequency sampling 

In practice, we define 

V = vi and 



I [0,JV/2-l] 



where Tobs is the observing time, N the number of equally 
spaced sampling points in the data set and I the bin 
number in the Fourier spectrum 



3.1 



we 

th 



Following Sec 

define Cn(^) ^ the discrete Fourier response of the n 
harmonic of a signal of constant amplitude emitted by a 
pulsar whose phase is modulated by circular orbital mo- 
tion. We also define its power spectrum P|^^„(£) = |Cn(^)P- 
Henceforth for convenience, various quantities such as or- 
bital frequency i^orb will be expressed in units of Fourier 



bins instead of frequencies. Thus, Vorh can also be referred 
to as Np^^^ = t'orb Tobs in the context of a particular ob- 
servation of length Tobs- The previous results then imply 
that 



1, 



Pi^,„ is a spiky A^p^^,^ -periodic sequence with a modu- 
lated amplitude. The periodicity is now expressed in 
units of Fourier bins. Also note that 



(11) 



where Np^^^ ~ Vpsr Tobs is the number of pulse periods 
in the observation. 
2. The Fourier response C„ is compact on the interval 
/(</?„, A^p^^j^) about the n^^ harmonic. It corresponds 
to the Doppler shift interval of the n^^ harmonic of the 
pulsar signal expressed in units of Fourier bins and can 
be approximated by 2 Np^^,,_^ ipn (see Fig. ^) . 

A A^-point Fourier Transform is elegant when the pro- 
cessing scheme is cast as a spectral decomposition in an 
iV-dimensional orthogonal vector space. However, such re- 
markable efficiency is compromised in practice because the 
number of orbital periods and pulse periods contained in 
an observation are not integers. 

3.2.2. Spectral leakage 

From the continuum of possible frequencies only those 
which coincide with the basis of discrete frequencies vi 
will project onto a single complex vector in the Fourier 
spectrum. All other frequencies will exhibit non zero pro- 
jections on the entire basis set. This is often referred 
to as spectral leakage and results from processing finite- 
duration signals whose frequencies are not periodic in the 
observation window W . This leakage can strongly degrade 
a signal's Fourier response. 

Multiplicative weighting windows can be used to minimize 
such an effect (Harris 1978). Such windows proved very 
useful for multiple-tone signal detection. However, they 
convey other intrinsic properties that make the detection 
of noisy signals difficult: 

— The Equivalent Noise Bandwidth^ increases. 

— The independent noise samples become correlated and 
the estimation of the probability density function of 
the noise consequently increases in complexity. 

— If the windows are applied to overlapping partitions 
of the data sequence, the noise samples become even 
more correlated. The work load also increases. 

Hence we compute the Fourier response given by the 
straight FFT of the data set. We can now rewrite Eq. 



^ The Equivalent Noise Bandwidth is the width of an ideal 
rectangular filter which would accumulate the same noise 
power from white noise as the window function's kernel with 
the same peak power gain. 
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(|) for the case of a discrete Fourier Transform as s„ {£) 
An{e) * [C«W * Si£ -nNp^J] where 



Int[-E(v)] 
fc=-Int[_E(ip)] 



(12) 



C„ can be calculated for virtually any window W pro- 
vided its Fourier response is known. In that sense, Eq. 

represents the analytical solution to the problem of 
Fourier analysis of noiseless binary pulsar signals in the 
solar system barycentric frame. In practice, W is defined 
by a Boxcar function whose Fourier response is given by 



VFsW =robse-^"^sinc(^) , sine (£) 



3.2.3. Computational cost 



sin {tt£) 
n£ 



(13) 



Obviously, C,n can be computed using Eq. (|l^). The sine, 
cosine and Bessel functions can be pre-computed and 
stored in static lookup tables. The computational cost of 
estimating ^„ over its compact interval for the simple case 
of the Boxcar function amounts to ~ 50 Np^^^^E{ipn)'^ op- 
erations. 

Alternatively, C„ can be evaluated numerically by comput- 
ing a Fast Fourier Transform (FFT) of a simulated data 
set of length K = Int [log2 li^Pn, ^Po,b) ] defined by 



Cn(fc) 



exp 



2'Kjn 



K 



K 



j(pn COS 2ttNp^ 



K 



(14) 



where fee [0 ; X — 1] . The spectral leakage factor is defined 
by 77 = Min[Frac[z^psi.T'obs], 1 - Frac[i'psrrobs]] such that 
77 e [—0.5 ; 0.5]. The number of operations required in that 
case to estimate Cn reduces to ^ 2 if (5 -I- logjif). We 
clearly favor the latter solution as the computational cost 
is much smaller. 

3.2.4. Additive noise contribution 

Let us now consider the effect of including additive white 
noise with the pulsar signal. The variance of the noise is 
assumed equal at all frequencies. Since the noise is additive 
and the Fourier Transform operator is linear, the noisy 
power spectrum is given by 



(15) 



where x{£) represents the Fourier decomposition of the 
noise samples at the frequency corresponding to bin £. 
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Fig. 3. Efficiency factor 7a obtained from numerical sim- 
ulations. The plot shows 5000 binary pulsar systems with 
random parameters. The distributions are uniform in the 
ranges given by Np^^^ e [2.0; 20.0], G [1.2; 9000], 
rye [—0.5; 0.5] and (j)orb & [— 7r;7r]. 

A{£) denotes the cross-terms between the signal and noise 
complex amplitudes : 



A = 2 [n{Sn) X 3?(i) + 3(S„) X 3(i)] 



(16) 



For a wide range of probability distributions, real (5R) 
and imaginary (3) parts of the Fourier decomposition 
of the noise samples are gaussian distributed (Jenkins & 
Watts 1968), hence \S:{£)\'^ will, with proper normaliza- 
tion, be chi-squared distributed with 2 degrees of freedom. 
The proper normalization factor P for a Discrete Fourier 
Transform (DFT) of length K defined by 



-2-iTjik/K 



(17) 



fc=0 



depends on the initial noise probability density func- 
tion. This function is gaussian in radio observations 
and poissonian in X-ray observations hence we find 
/3 — yjK/2. When the noise is chi-squared distributed, 
we find (3 = V2K. 

Because of the orthogonality of the Fourier decomposition, 
3? and 3 are independent. Therefore, A(^) is gaussian dis- 
tributed with mean fi{e) = and a{£) = ^/tp^^jJ). The 
statistical expectation of A consequently reduces to zero. 

3.3. Summary 

We have now completed our Fourier analysis of a peri- 
odic signal whose phase is modulated by circular orbital 
motion. For an observation of finite duration, we can sum- 
marize our results as follows : 



The basic analytical solution is given by Eq. (h2 
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— The power spectrum is iVp^^^j,^ -periodic over a finite fre- 
quency range about each harmonic. 

— The phase in Fourier space is unpredictable a priori 
because Np^^^^ and (path are unknown. 

— A numerical simulation including an FFT is more effi- 
cient than a direct computation of Eq. (jl^) . 

~ The statistical properties of white noise do not favor 
any particular periodicity. The particular signature of 
a pulsar signal whose phase is modulated by a circular 
orbital motion consequently provides an opportunity 
for an incoherent detection of weak pulsar signals in 
broadband noise. 

We have considered pulsars in circular orbits. The case 
of eccentric orbits is extremely difficult to solve analyti- 
cally. The envelope of the Fourier response becomes very 
complicated and can not be described by simple Bessel 
functions. However, the sideband spacing still remains the 
same and therefore produces a periodic sequence as in the 
case of circular orbits. 

We have not considered non-uniform sampling as may 
happen for long observations due to Doppler correction 
for the observatory, either earthbound or satellite-borne. 
We assumed that the time series has been resampled with 
equal intervals for an observation in the solar system 
barycentric frame. 

In Sec. ^, we briefly present three pulsar search tech- 
niques widely used in radio/X-ray astronomy and in- 
vestigate their limiting sensitivity to binary pulsars in 
close circular orbits. In Sec. ^, we use the results of the 
present section to develop a very efficient implementa- 
tion of the Phase Modulation Searching method (Ransom 
1999). Eventually, we compare the efficiency of these re- 
covery techniques and discuss the detectability of a range 
of binary pulsars in ultra-short circular orbits. 

4. Searching for pulsars in close binary systems 

When a pulsar signal is present in a time series, which in 
the case of radio observations has already been corrected 
for the dispersion due to its passage through the interstel- 
lar medium, it manifests itself as a set of discrete harmon- 
ics in the Fourier spectrum. The number of significant har- 
monics depends on the duty cycle of the pulsar (fractional 
on-time of the pulsar signal). Whereas ordinary pulsars 
have a typical duty cycle of about 4% (indicating about 
12 significant harmonics in the spectrum), the millisecond 
pulsars seem to have a duty cycle of about 10% (about 5 
harmonics). In the population of X-ray pulsars, the High 
Mass X-ray Binary pulsars have complex pulse profiles 
and thus have a few significant harmonics. However, some 
Low Mass X-ray Binary pulsars, including the only known 
millisecond pulsar among them, show nearly sinusoidal 
pulsations. Naturally, whenever we have more than one 
harmonic, "harmonic-folding" (Bhattacharya 1998) will 
improve the probability of detection when compared to 
looking for statistically significant isolated harmonics in 
the power spectrum. 



In the case of solitary pulsar signals, harmonics are es- 
sentially confined to one bin whose power is decreased by 
sinc^(Min[|Frac[n77] 1, 1 — |Frac[7i77] |]), where n is the har- 
monic number ; it is important to account for the spectral 
leakage factor 77 because we will compare different recov- 
ery techniques that are not equally affected by such an 
effect. In a way similar to that of Johnston & Kulkarni 
(1991), we define the efficiency factor 7, such that a bi- 
nary pulsar must, on average, carry 7^^ more pulsed fiux 
density than its virtual solitary counterpart to be detected 
with the same significance level. We now present several 
methods and calculate their 7 either from theoretical pre- 
dictions or numerical simulations. We will compare and 
discuss their performances in fuller details in Sec. 0. 



Method A 

Consider a signal emitted by a pulsar in a binary sys- 
tem with a circular orbit. This signal is corrected to 
the solar system barycenter if necessary and uniformly 
resampled. Then, we simply search for statistically sig- 
nificant peaks in the discrete power spectrum - hence- 



forth referred to as P-space. As seen in Sec. 3.2.1 



a given harmonic number n will be spread out over 
I{LPn, Np^^-^^) bins. The efficiency factor 7a is plotted 
as a function of ip„ for 5000 systems with random pa- 
rameters in Fig. 1^. The distributions are uniform in the 
ranges given by Np_^^^^ € [2.0; 20.0], ipn e [1.2; 9000], 



7/ e [-0.5; 0.5] and 



rb € 



-7r;7r]. Clearly, 7a is 



strongly correlated with (p„ . The dispersion about the 
least-square-fit line 



7A ^ [0.74 ±0.1 



(18) 



is mostly caused by spectral leakage effects. Therefore, 
7a is maximized when Np^^^^ and 77 are close to integer 
values. 

As ifn oc ai and 7a varies little with orbital period, we 
can see in Fig. ^ that close binaries are easier to detect 
than wide binaries. This is because a pulsar with a 
relatively small (pn experiences less acceleration. 
In the worst case, sensitivity limits are decreased by a 
factor of about 30 when compared to a solitary pulsar. 

Method B 

Numerous authors have developed techniques to cor- 
rect for the Doppler smearing caused by the orbital 
motion. However, a complete coherent recovery in- 
volves searching in a parameter space defined by a 
large range of keplerian parameters. This is com- 
putationally unfeasible for the present day technol- 
ogy. As a compromise, some authors (Middleditch 
& Priedhorsky 1986; Anderson et al. 1990; Wood et 
al. 1991; Camilo et al. 2000; Johnston & Kulkarni 
1991) have attempted a constant acceleration correc- 
tion. Our simulations showed, in good agreement with 
Camilo etal. (2000), that almost 100% of the pulse 




Fig. 4. The binary system is described by Vpsi — 1 KHz, Poib = 1-5 hours, M2 = 0.8 Mq and (/)o, b — 53° so that 
ipi — 4571 radians, Np^^^^ — 5.34 and 77 ~ 0.23. Upper plots show raw power spectra as obtained from a direct DFT 
of a noiseless/noisy observation (left/right respectively). Lower plots show Q. The pulsar is detected at as 
shown by feature (1). The second harmonic is indicated by feature (2). Features (la) and (2a) are essentially caused 
by beating of Np^^^^ with bins in 7^-space. The relative position of (la)/(l) and (2a)/(2) is given by l+Frac[ iVp^^^ ]. 
Multiplicative weighting windows could almost remove such features caused by spectral leakage effects (see Sec. |3.2.2 ). 



strength is on average recoverable for all orbital phases 



if Nf 



1/7. The efficiency of constant acceleration 



searches can therefore be expressed as 



7b 



(19) 



Clearly, constant acceleration searches perform well 
when the observing time is much less than the orbital 
period. The efhciency of such a method can be ex- 
pressed in various ways. We have estimated 7b in order 
to maximize the recovery for the smearing due to the 
time-dependent Doppler effect experienced by the bi- 
nary pulsar. This method basically allows for searching 
the whole time series with negative and positive accel- 
eration values, including zero acceleration. Therefore 
in practice, it includes Method A which may itself 
prove more sensitive when A^p^^^ large. However, 
Method B and 7b will henceforth refer to Eq. (|l9|) 
as aiming for a "non zero" constant acceleration cor- 
rection. 

Method C 

Another method widely used in X-ray astronomy con- 
sists of dividing a long time series into sub-segments 
in order to reduce the frequency resolution so that the 
Doppler shifted pulsar signal remains essentially con- 
fined to one bin in T'-space (van der Klis 1989). The 
A^s short power spectra are then incoherently summed. 
The stacked power spectrum is subsequently searched 
for statistically significant peaks. 

Using confidence levels instead of signal-to-noise ratio 
(SNR) is of great importance here because the noise 



distribution itself is modified in this process. The effi- 
ciency factor 7c can therefore be measured as follows 

D f 

: let us define Xa as the inverse cumulative density 
function for a chi-squared distribution with D f degrees 
of freedom at the confidence value a {Df = 2 Ag). By 



analogy with Sec. 3.2.4, we can define X = l^fc 



so that Prob(X > Xa^) — ^- The average expectation 
of X is given by Xxi2- When Df is large, the noise 

distribution approaches a gaussian and xf/2 coincides 
with the average value of X. Also, when Df is large, 
the noise can be considered additive in power instead 
of in Fourier amplitude (Vaughan et al. 1994; Groth 
1975). Hence, 7c is given by 



7c 





- xl/2 


2JV3 


-^1/2 



(20) 



7c must be evaluated at discrete values given by A's — 
2P,p G N. In a sense, such an approach is very simple 
as 7c depends only on I{ipniNp^^^^), but we can not 
immediately recover any orbital parameter. However, 
this method can be very powerful (see Sec. || for a fuller 
discussion). 

5. Partial coherence 
recovery technique 

We now describe our Partial Coherence Recovery 
Technique. We start with an observation of radio/X-ray 
pulsations which has been Fourier Transformed as in 
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Table 1. Fourier spectrum statistics. The kurtosis is mea- 
sured on real and imaginary Fourier samples obtained 
from an exponentially distributed noise. A normal distri- 
bution has Kurtosis = 3. Higher values indicate leptokur- 
tic distributions. 



10" 



K 



16 



32 



128 256 512 1024 



Kurtosis 3.83 3.51 3.25 3.10 3.05 3.02 3.01 



Method A. Somewhere in T'-space lies a set of A^p„,b- 
periodic sequences centered on harmonic frequencies of 
Np This periodicity, often buried beneath the noise, 
can be detected by computing a short DFT of Vn of 
length K ^ I(ipn,Np^^^) bins. Subsequently, we define 
S = l^nPj hereby defining Q-space as some sort of 'in- 
complete' T-space (T-space is obtained by an inverse DFT 
of Fourier amplitudes whereas we perform a forward DFT 
of Fourier powers). Searching for statistically significant 
peaks in Q would reveal the presence of a sufficiently 
strong pulsar signal at if Np^^,^^ > 2 (see Shannon's 

Sampling Theorem - Shannon 1948). Thus, a detection is 
made as shown in Fig. ^ Both pulsar and orbital phases 
are lost in Q-space, hence we call this technique Partial 
Coherence Recovery Technique (henceforth referred to as 
PCRT and Method D). 

The following sections are devoted to pulsar signal detec- 
tion and estimation, where detection is the task of de- 
termining if a pulsar signal is present in an observation, 
while estimation is the task of obtaining the values of the 
parameters describing the signal. 

5.1. Detection 

Since Vn is localized in 'P-space and its extent unknown a 
priori, we must compute many DFTs at all possible har- 
monics of Nppsr with various widths K. We subsequently 
use a set of sliding windows W and define Wk.i as the 
window of size K centered on bin £. A sufficiently strong 
pulsar signal will be detected when Wk/ contains Vn, i.e. 
Np^^, i,ndK>Iiipn,Np^J. 



10" 



10' 



10 



10" 



10 



Fig. 5. The dot-dashed curve (Method A) is given by Eq. 
(|lf). The solid curve represents Method D (PCRT). 

up more noise samples than we would expect on the ba- 
sis of a gaussian distribution. The lower limit to assign 
to K is however somewhat arbitrary as one has to define 
a criterion of exclusion where an estimator deviates from 
gaussian statistics; we chose Kynin — 2*. 
The upper limit is determined by Max[/((^„, iVp^_.j^)]. It 
corresponds to "the" extreme astrophysical case where the 
number of bins over which the signal power is distributed 
is maximum. Let us consider a hypothetical 0.5 ms pulsar 
of 1.4 Mq with a companion star of 1.4 Mq in a binary 
system with a circular orbit of 3.728 hours; i = 90 degrees, 
-^forb = 2 and n = 2. We use the second harmonic as we 
consider an extreme case and therefore aim to maximize 
(fin. The orbital velocity of the pulsar Wmax — 279 km/s and 
the second harmonic power spreads over about 200 x 10'' 
(<2^*) bins. We therefore take an upper limit i^^max of 
218. 

The number of values possibly taken by K is given by 
A^A' l + log2 J^max — log2 -K'min- In the case we considered 
above, = 11. 



5.1.1. Sizes of the short DFTs 

Since Pc^^n is compact, K must be ~ I{'.pn,Np^^^) bins to 
achieve optimal detectability. Computing DFTs of random 
sizes is not achieved efficiently with modern algorithms, 
so we must compromise our sensitivity to minimize the 
computational effort by computing the DFTs using FFTs 
of dyadic length if = 2P,p e N. 

We need a well-behaved noise probability density function 
to compute appropriate detection levels. The noise distri- 
bution in Q-space will approach gaussian statistics when 
K is large enough (Central Limit Theorem). In Table |l|, we 
show that such a distribution is leptokurtic for small val- 
ues of K. This implies that more probability is distributed 
in the tail, so a larger number of high noise values must 
be expected. In other words, a search code would pick 



5.1.2. Sliding windows 

The ideal situation occurs when the power spectrum of 
'WK,t can be computed at all possible frequency bins £ 
where a pulsar signal is possibly expected. Let us denote 
that number by L. If the sampling time is 0.1 ms and 
minimum pulsar period 0.5 ms, L « N/h. The number of 
operations required to compute all FFTs is given by (see 
Appendix C.l) 



OfFT ^ 2^^^LX^i„ [l0g2 i^min + N K - 2] 



(21) 



In our example where Nk — 11 and i^min = 256, Offt ~ 
9i X 10^. If A = 228, Offt ~ 4.8 x lO^"* which requires 
about 5.6 days of computations on a 1 GigaFLOPSQ ma- 

^ 10^ floating-point operations per second (FLOPS). 
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Fig. 6. The dot-dashed curve (Method B) is given by Eq. 
( p^ ) . Note that 7b scales with A^p^,^ • Solid curves (Method 
D) correspond to K = 2^ . . .2^^ (upper curve is = 2^). 
Higher 7 implies better sensitivity. Method D is superior 
to method B in nearly all cases. 



chine. This is prohibitive. An alternative is to decrease 
the fraction of overlap g between consecutive Wk/ where 
g=l~p/K,pG[l: i^min — 1] (note that g is a. frac- 
tion, whereas p denotes a number of bins) . The number of 
operations then reduces to (see Appendix |C.2[) 



O 



L 



FFT 



1-g 



Nk l0g2 i^min + 



Nk{Nk-1) 



(22) 



li g = 0.87, Offt ~ 5.9 x 10^° which requires about 1 
minute of computations on a 1 GigaFLOPS machine. 

A sliding window of length K must catch Pf,n with 
I{(pn, Np^^^^) < gK. For such signals where the proba- 
bility of being entirely swept up by Wk.^ is 1, at least 
1-g fraction of the sliding window will convey only noise. 
In that sense and because Wk/ is evaluated for a small 
number of sizes Nk , we can say that our algorithm suffers 
from our limited computation power. 
The number of correlated spectral samples eventually ob- 
tained after processing L bins with a single sliding win- 
dow of length K amounts to L/2{1 — g). Overlapping 
windows cause another intrinsic problem, namely the cor- 
relation of random components. If Wk,i is defined with 
a Boxcar function, successive VVif/ will be 100 x g % 
correlated. This implies that the number of uncorrelated 
spectral samples obtained after processing L bins with a 
single window of size K drops back to L/2. 

5.1.3. Detection levels 

y^K,i contains a chi-squared distributed broadband noise, 
and may also include a jVp^^ -periodic sequence caused by 
a signal. As seen in Sec. 3.2, such a sequence is very spiky. 



Table 2. Harmonic folding ambiguity. L/Nh is given by 
Eq. ( p3| ) and TV = 2^^. Na denotes a number of combina- 
tions. 



H 



6 



7 



8 



9 



L/N„ 
Na 



12 
3 



24 
5 



40 
9 



60 
13 



84 
17 



112 
23 



144 
31 



180 
37 



This comb of sidebands essentially decomposes into H ~ 
Int[iVp^^j^/2] harmonics in Q-space if A'^p^^,^ > 2 (Shannon 
1948). The presence of a sufficiently strong binary pulsar 
signal would be detected in Q at harmonics of its orbital 
period as shown in Fig. IJ. 

When Np^_^^ > 4, it becomes useful to fold successive 
harmonics in Q to improve detection limits. Harmonic 
folding (Bhattacharya 1998) denotes the process of inco- 
herently gathering harmonically related signal peaks in a 
power spectrum. For instance in Fig. ^, we clearly see that 
adding the two harmonics, namely features (1) and (2), 
will greatly improve the significance level of the detection. 
Folding means incoherent addition of H independent chi- 
squared distributed samples resulting in chi-squared dis- 
tributed noise with 2Ji degrees of freedom. Since a de- 
tection has a statistical significance for a particular noise 
distribution, the number of uncorrelated spectral samples 
N}i to be searched for a pulsar signal must be dealt with 
separately as a function of . In other words, the signif- 
icance of a power obtained after folding 2 or 3 harmonics 
is different because the noise distribution itself depends 
on the harmonic fold number. Therefore, we define Nn as 
the number of samples that can be folded a maximum of 
H times in Q-space {H = 1 means no folding since we 
consider only harmonic number 1). A simple dichotomy 



gives 



1 

li 



H+1 



L 



2H{H +1) 



(23) 



where, for a given sliding window of size K, L/2 represents 
the total number of uncorrelated powers in Q-space and 
the term within square brackets denotes the fractionnal 
amount of samples that can be folded a maximum of H 
times. 

Another complication arises because each frequency chan- 
nel has a finite bandwidth. We can incoherently add bin 
number 10 with bin number 20, but also with bin num- 
ber 19 or 21 depending whether ?/ < —0.25, jryj < 0.25 or 
1] > 0.25. This ambiguity becomes extremely complicated 
when H increases and requires a large amount of book- 
keeping. We define Na as a measure of such an ambiguity. 
In other words, Na represents the number of possible com- 
binations in folding H harmonics. As shown in Table ^ 
Na increases with harmonic number H. A comprehensive 
search must cover all possibilities. Therefore, the signifi- 
cance level of a detection obtained after folding must take 
into account the number of combinations Na as well as 
the appropriate noise distribution. 
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The detection level at harmonic number H is consequently -i o° 
given by X^"naNhH^,^^ where i?max = 9 in our case. We 
used a confidence level a of 1% to estimate 7d- 



5.2. Estimation 

Once a significant peak is detected and locafized both in 
■P-space and Q-space, we can perform efficient complex 
cross-correlations of length ~ /((^„, A'^p^^.^) with a rela- 
tively small number of functions C,n given by Eq. (p^. 
Such a coherent recovery of a suspected pulsar signal is 
very computationally efficient and extremely sensitive. As 
we showed in Sec. 3.2.3, Cn is defined by four parameters 



1. Porb is obtained directly from Q where a pulsar signal 
has first been suspected as in Fig. |^. The accuracy to 
which this parameter is known improves with (p„ and 

2. ipn can be approximately estimated given Np^^^^ and 
the length K of the current sliding window. However, 
a more accurate estimation can be obtained efficiently 
by performing sliding FFTs of various lengths about K 
in order to maximize the SNR at the frequency corre- 
sponding to 1 /Np^^^^ . Since we assumed circular orbits 
Pi,n is symmetric, so this process also yields an ap- 
proximate value of Vpsr • 

3. (^orb can be obtained by performing complex cross- 
correlations in T-space. Since we work on a small 
number of Fourier amplitudes (as I{(pn, Np^^^^) ^ 
2Np^^^^^ ip„ bins), we can inverse Fourier Transform Q 
and complex cross-correlate with our reference func- 
tion C,n given by Eq. (|lj) and (poih = 0. The initial 
orbital phase 0orb can be easily recovered if Np^^,^^ is 
known. In practice z/psr is approximated when comput- 
ing Cn- 

4. 77 gives fpsr within about half a frequency bin of pre- 
cision. This parameter can be obtained by complex 
cross-correlating Q with ^„ using previous estimates 

of Afp„b' 'fin and (/)orb- 

The pulsar frequency and all measurable binary param- 
eters are therefore known at the end of this procedure. 
Such a recipe holds for pulsars in circular orbits. If the 
periodic sequence in P-space is caused by a pulsar in ec- 
centric orbit, another procedure must be used and further 
processing is required. We have not investigated such a 
case as our aim is to give a complete recipe for searching 
for pulsars with circular orbital motion. The case of ec- 
centric orbits should however receive full attention in the 
context of the Phase Modulation Searching method. 

6. Discussion 

We have considered binary pulsars in circular orbits in Sec. 

In Sec. ^, we derived the analytical solution to the prob- 
lem of Fourier analysis of such pulsar signals in the solar 
system barycentric frame, in the continuous and discrete 



10" 



10^ 



10' 
K 



lo-' 



Fig. 7. The dot-dashed curve (Method C) is given by Eq. 
( pO| ) . 7c depends only on K 2 Np^^^ (fn . Solid curves 
(Method D) correspond to TVp^^j, G [2H;2{H + 1)[ where 
_ff = 1 . . . 9 is the harmonic fold number defined in Sec. 
|3.1.3 , The upper curve corresponds to iJ = 9. The 7 values 
are estimated at K — 2^, p £ N. 



regimes. We found that each harmonic in the power spec- 
trum of such signals spreads over about 2 Np^^^^ ipn bins, 
within which it is Np^^^^ periodic. We suggested two ways 
of modelling the Fourier response of such pulsar signals 
and estimated the computational cost of these methods. 
In Sec. ^, we have estimated the detectability of such pul- 
sars (assuming that all parameters were unknown) using 
three widely used methods (Method A : standard Fourier 
analysis - Method B : constant acceleration correction - 
Method C : stack searches) for which we derived efficiency 
factors 7. In Sec. ||, we presented a detailed description 
of our Partial Coherence Recovery Technique (Method D) 
which takes advantage of the periodicity in the power spec- 
trum of a binary pulsar signal (see Sec. |^) to implement 
a pulsar search algorithm based on the Phase Modulation 
Search method (Ransom 1999). In this section, we present 
and discuss the sensitivity of all four techniques to close 
binary pulsars of different types with various companion 
masses and orbital periods. 

We have simulated 2^^^ random binary systems with 

ipi e [1; 10000] 

A^P„., e [1.5; 20] 

(port e [-7i';7r] 

?7 e [-0.5; 0.5] 

and calculated 70 for each system individually. We com- 
pare our results with Method A, B and C on Fig. ||, || and 
0. We recall that the sensitivity of pulsar search algorithms 
usually expressed in units of milli-Jansky is inversely pro- 
portional to 7. 

As shown in Fig. ||, 7d/7a^ 3 for all ipn and both meth- 
ods are strongly correlated with However, 70 varies 
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Fig. 8. Contour map of (pi as a function of orbital period 
Porb and companion mass M2. We assumed i = 90° and 
z/psr = 250 Hz. We recall that ipi (x Vpsr sin(i). 



relatively little with Np^^^ as we can see in Fig. |g. Fig. 
^ and 1^ compare state-of-the-art methods in recovering 
smeared pulsar signals. Method D is significantly better 
than Method B when searching for very close binary sys- 
tems with ultra-short orbital periods for which Np^^^^ is 
large (e.g. Np^^^^ = 36 for Porb = 10 minutes and Tobs = 6 
hours) and relatively small ifn- Furthermore, 7b in Fig. 
^ represents an ideal case where it is possible to run a 
constant acceleration code on about one seventh of Porb , 
which is rarely the case as the orbital period is usually 
unknown. Therefore, we expect 7d/7b to be even larger. 
Another important consideration concerns the processing 
time. In Method B many long FFTs are taken repeatedly 
for each positive and negative acceleration trial for vari- 
ous data lengths, whereas in Method D only one long FFT 
is performed. The gain in processing time is considerable, 
especially when the range of orbital accelerations is large. 
Method D (PCRT) does not perform well when Np^^^ < 2 
because a detection then relies on an aliased frequency 
component (Shannon 1948). 

Fig. 1^ demonstrates that Method D is always more sensi- 
tive than Method C. The reason is as follows : as seen 
in Sec. ^, a harmonic manifests itself in P-space as a 
comb of sidebands separated by Np^^-^^ — 1 bins. These bins 
carry mostly noise and are co-added de facto in Method 
C. Whereas in Method D, we perform a spectral analysis 
which basically ignores them. However, Method C must 
be regarded as a useful method when compared to other 
standard searching techniques. Its good performance re- 
sides in the distribution of the noise itself, which becomes 
gaussian because many values are co-added, while, due to 
the decrease in frequency resolution, one deals with many 
fewer points hence improving the significance level of a 
detection. Such a method is also computationally very ef- 
ficient. A clear drawback however resides in its lack of 



information about the binary system and poor frequency 
resolution. 

In order to estimate the detectability of pulsars in ultra- 
short binary systems, we have plotted tpi as a function of 
Porb and companion mass M2 for a 4 ms pulsar in Fig. ^ 
This plot can be used for determining ipi for any particular 
pulsar. For example, let us consider a 8 ms pulsar in a 
binary system with Porb = 1 hour and M2 = 0.1 M©. The 
value of (pi is 128 * 4/8 = 64 radians, where 4 corresponds 
to the reference pulsar period in Fig. || and 8 corresponds 
to the pulsar period we consider in milliseconds. Of course, 
ipi is a maximum value since we assumed an inclination 
angle of 90 degrees but we recall that ipn scales directly 
with sin(i). Let us assume that Tobs = 6 hours, hence 
-^foib = 6. Looking up the 7 values for Method A, B, C 
and D in Fig. I, I and @, we see that this particular binary 
system would be best detected using Method D (PCRT) as 
7d > 7c > 7A > 7b (see Sec. || for an accurate definition 
of 7A and 73). 

Of particular interest is the detection of pulsars with large 
companion masses such as Neutron Star - Neutron Star 
and Neutron Star - Black Hole systems in circular orbits. 
In such systems, millisecond pulsar signals are smeared 
over many bins in the Fourier spectrum because ifn is 
large. These systems are characterized by large accelera- 
tions which are usually not searched for in constant accel- 
eration searches due to the large increase in the number of 
acceleration trials, which tremendously increases the pro- 
cessing time. Stack searches could be very useful in such 
a context. However, the PCRT is, on average, the most 
efficient and sensitive method available for detecting such 
extreme systems, with the proviso that the orbital period 
is short so that -/Vp^,.,^ is greater than 2. 
By direct and detailed comparison with other methods 
for correcting for the period variation caused by orbital 
acceleration, we find the PCRT is a very efficient method 
for detecting fast rotating radio/X-ray pulsars in binary 
systems, especially when A^p^^^ is large. However, it re- 
quires long observations, and therefore is not well suited 
for direct application to all-sky pulsar surveys where in- 
tegration times are usually less than half an hour. On the 
other hand, the method is ideally suited to long observa- 
tions of globular and open clusters, steep spectrum point 
sources and X-ray point sources. 
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Appendix A: 

Main symbols defined in the text, with description and section of first use. 



Symbol 



Description 



Section 



7 
H 

K 
L 

V 
N 

Na 
Nh 
Nk 

Pn 
Q 

Q 
P 

Cn 

w 

J^orb , Np. 



Inverse cumulative density function for a chi-squared 
distribution with Df degrees of freedom. 

Efficiency of a binary pulsar search method ( 0.0 < 7 < 1.0 ). 

Harmonic fold number in Q-space. 

Width of in bins (about 2 A'^p^,,, bins) . 

Length of a sliding window (Wk) in bins. 

Number of spectral samples searched for a pulsar signal in P-space. 

Spectral leakage factor : rj £ [—0.5; 0.5] . 

Number of equally spaced sampling points in T-space. 

Number of possible harmonic fold combinations. 

Number of samples which can be folded exactly H times. 

Number of sliding windows Wk used in the analysis. 

Number of sub-segments used in Method C. 

Power spectrum of (~ liVn, Np^^^^) bins wide). 

Noisy Pf,„. 

Power spectrum of Vn- 

Percentage of overlap between consecutive Wx. 
Number of overlapping bins between consecutive Wk- 

Discrete Fourier response of the n**^ harmonic of a noiseless binary pulsar signal 
of constant amplitude, shifted back to the origin of the frequency scale. 

Observing window e.g. Boxcar (Wb). 

Sliding window of length K. 

Half the phase rotation experienced by the harmonic of the pulsar signal 
during the orbit, in radians : = limVp^,- [ai sin(i)/c]. 

Initial orbital phase in radians, measured from superior conjunction of the pulsar. 
Observed pulsar rotation frequency in Hz and Fourier bins, respectively. 
Observed pulsar orbital frequency in Hz and Fourier bins, respectively. 



3.2.1 


3.1 




5.1 





Appendix B: 

We derive the mathematical proof of Eq. (^. We shall demonstrate that 



00 

, jz COS " 



= i''Jk{z)e^^ (B.l) 



fc— — 00 

By definition, 

gjzcose ^ cos(zcos6i) + jsin(zcos6l) (B.2) 
where, according to Eqs. (9.1.44) and (9.1.45) in Abramowitz & Stegun 1974, 



00 

cos(zcos6') = Jo(z) + 2 Y^{-IY J2p{z) cos{2pd) 
p=i 



(B.3) 
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sin(zcos6') = 2 ^(-If J2p+i(2) cos( [2p + 1] 0) . 

p=0 



Since j = (—1)^/^, Eq. (B.4) can be rewritten as 



(B.4) 



sin(zcos6') = - ^ {-ly J2q{z) ,cos(2q6') 

^ 9=1/2 



where q = p + -k- Combining Eqs. (B.3) and (|B.5|) in (B.2) yields 



(B.5) 



gj.cose ^ Jo(^) + 2 cos(fc( 



k=l 



Since cos{ke) = [e^'^'^ + e-^''^] /2, we get 



k=l 



According to Eq. (9.1.5) in Abramowitz & Stegun 1974 : Jfe(z) = j J^k{z). Therefore 



k=l 



which simplifies to 



k— — oo 

We note that in Eq. j'' has been replaced by its equivalent exponential notation e^^^/'^ . 

Appendix C: 

C.l. Determination of Eq. 



(B.6) 



(B.7) 



(B.J 



(B.9) 



Let us consider a given frequency channel and compute Nk FFTs of length K S [ivTmin ■ • ■ -f^max] ■ The number of 
operations required for an EFT of length K is given by i<'log2 K . Therefore, 



OfFT = i^minl0g2 i^min + ■ ■ ■ + [2'''<-^K^^] log^ [2^'<-'Kn 



(C.l) 



which can be rewritten as 



Cfft — Kn 



Nk-I Nk-1 



log2(i^.„i„) ^ 2^+ ^ 2'^-fc 



fc=0 fc=0 



The summations can be explicitly computed so that 



OfFT = i^min [ log2(i^mi„) [ 2^^' - 1 ] + [ 2^^" (A^^ - 2) + 2 ] ] 



Assuming 2 ^ 1, Eq. (C^) subsequently reduces to 



OfFT ^ 2^^Kmin [l0g2(ifmin) + Nk - 2] 

Naturally, Offt niust be multiplied by L hence Eq. (pT]). 



(C.2) 



(C.3) 



(C.4) 
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C.2. Determination of Eq. (22) 

We now consider p > 1. The number of frequency channels at which a FFT of length K is evaluated reduces to 
L/K{1 - g). Therefore, 



Cfft — 



This equation can be rewritten in a simple form as 



2^--l/fmin(l - g) 



(C.5) 



O 



L 



FFT 



1-g 



Nk-1 
k=0 



(C.6) 



Subsequently, we obtain Eq. 



O 



FFT = 



L 



Nk log2 



Nk{Nk - 1) 



(C.7) 
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